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Recent studies show that spherical motile micro-organisms in turbulence subject to gravitational 
torques gather in down-welling regions of the turbulent flow. By analysing a statistical model we 
analytically compute how shape affects the dynamics, preferential sampling, and small-scale spatial 
clustering. We find that oblong organisms may spend more time in up-welling regions of the flow, and 
that all organisms are biased to regions of positive fluid-velocity gradients in the upward direction. 

We analyse small-scale spatial clustering and find that oblong particles may either cluster more or 
less than spherical ones, depending on the strength of the gravitational torques. 


PACS numbers: 05.40.-a,47.63.Gd,47.27.-i,92.20.jf 

Patchiness in suspensions of micro-organisms is fre¬ 
quently observed on a range of spatial scales. The under¬ 
lying mechanisms differ, depending on the properties of 
the micro-organisms, and upon the spatial scale. Patch¬ 
iness can be caused by density stratification and ver¬ 
tical shears [1], by predator-prey cycles, or by interac¬ 
tions between the organisms and water-column gradients 
- in light, chemistry, turbulence, and in hydrostatic pres¬ 
sure [2]. Patchiness is important because many biological 
processes (mating, feeding, predation) rely on individual 
encounters [3], and the encounter rate is strongly influ¬ 
enced by small-scale number-density fluctuations. 

Gravitaxis may cause such inhomogeneities in the spa¬ 
tial distribution of motile micro-organisms. Density- or 
drag-asymmetries of the body give rise to torques affect¬ 
ing the swimming direction [4-6]. When the effects of 
gyrotactic torques and fluid-velocity gradients balance, 
inhomogeneities may form in the spatial distribution, as 
shown by the nricro-alga Chlamydomonas nivalis swim¬ 
ming up against a down-welling pipe flow. The rnicro- 
algae gather in the centre of the pipe where the down- 
welling velocity is largest [7]. Gyrotaxis may trap motile 
organisms in macroscopic shear gradients [8, 9], and fluc¬ 
tuating vorticity may cause patchiness [10]. This is con¬ 
firmed by recent direct numerical simulations (DNS) of 
motile, spherical micro-organisms in turbulence [11] re¬ 
vealing that the organisms are more likely to be found in 
down-welling regions of the turbulent flow, they ‘prefer¬ 
entially sample’ such regions. 

These results raise three fundamental questions that 
we address and answer in this Letter. First, how does 
shape affect the dynamics in turbulence of motile micro¬ 
organisms subject to gyrotaxis? In Ref. [11] the organ¬ 
isms were assumed to be spherical. Non-spherical or- 
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ganisms respond not only to turbulent vorticity but also 
to turbulent strain [12-15]. This causes passive rods 
to exhibit intricate orientational patterns on the sur¬ 
face of turbulent and other complex flows [16-18]. Also, 
shape strongly affects the trajectories of active particles 
in model flows [19-21], and recent DNS indicate that pro¬ 
late gyrotactic organisms cluster less than spherical ones 
when gyrotaxis is strong [22]. Second, where do the or¬ 
ganisms go in turbulence? Are there circumstances where 
the organisms may not gather in down-welling regions, 
or where other mechanisms of preferential sampling may 
apply? Third, the fact that orgamisms tend to gather in 
certain regions of the flow (preferential sampling) does 
not explain which mechanisms actually cause them to 
get in contact. To determine these one must follow the 
dynamics of two organisms that are initially very close 
together, and determine whether they tend to approach 
further or move apart. We refer to the resulting small- 
scale spatial fluctuations in the number density as ‘small- 
scale clustering’. 

Statistical model. To answer these questions we use a 
simplified model [7, 11, 22] for the translation and ro¬ 
tation of small axisymmetric active particles subject to 
turbulence and gyrotaxis: 

r = v = u(r,t) + v s n and h = u>(r, t) A n . (1) 

Dots denote derivatives w.r.t. time t, r is the particle 
position, and u is the flow velocity. Each particle swims 
with constant speed v s in the direction n of its symmetry 
axis (|ro = 1). The angular velocity of the particle is 

w(r, t) = (g A n)/(2SB) + Ul(r, t)+An A [S(r,f)n] . (2) 

The first term on the r.h.s describes gyrotaxis. The 
unit vector g points in the direction —e z of gravity, and 
SB is the reorientation time [7, 11]. It depends on the 
mass distribution within the particle, and on its shape 
through hydrodynamic resistance. The other terms on 
the r.h.s. of Eq. (2) represent the effect of the turbu¬ 
lent velocity gradients upon the particle orientation [12]: 
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fl = (V A u)/ 2, and S is the symmetric part of the ma¬ 
trix A of fluid-velocity gradients. The parameter A char¬ 
acterises particle shape: A = 0 for spheres, and A = 1 
for infinitely thin rods. Eq. (2) disregards turbulent ac¬ 
celerations. In most marine conditions this is an excel¬ 
lent approximation [23]. We model the dissipative range 
of turbulence by incompressible, homogeneous, isotropic 
Gaussian random functions with typical length- , time- 
, and speed-scales rj, r, uq [24]. This neglects inertial- 
range properties which may become important for parti¬ 
cles that are larger than the Kolmogorov length [25]. We 
note that the dissipative-range turbulent fluctuations are 
universal [26], but they are not Gaussian. We comment 
on this difference between turbulence and the statistical 
model below. 

There are four dimensionless parameters: the shape 
parameter A, the reorientation time = 38 /t, the 
swimming speed 4> = v s t/t), and the Kubo number 
Ku = Uqt/t]. We vary the parameters independently, 
keeping 38 constant as A is changed. Ku is a dimension¬ 
less measure of the correlation time of the flow. 

Our choice of the dimensionless parameters is dictated 
by the method (explained below). DNS employ differ¬ 
ent de-dimensionalisations [11]: 38 by the Kolmogorov 
time tk = l/y / Tr(AA T ) 00 , and v s by the corresponding 
Kolmogorov speed uk- Our dimensionless parameters 
translate to those used in the DNS as Tdns ~ Ku 'F and 
d?DNS ~ 4> / Ku. We expect that the statistical-model 
results become independent of Ku at large Ku and qual¬ 
itatively agree with DNS results [24]. 

Typical values of v s are given in Ref. [5]: 38 ~ 1 
5 s and v s ~ 0.1 -lmm/s. Typical ocean dissipation 
rates are e ~ l-10 2 mm 2 /s 3 for surface water [27], giv¬ 
ing Kolmogorov times, lengths, and speeds in the range 
tk ~ 0.1-1 s, r] ~ 0.3-1 mm, and Uk ~ l-3mm/s. These 
estimates yield Tons ~ 1 50 and 4 >dns ~ 0.03-1. In 
Ref. [28] smaller dissipation rates, e ~ 10 _4 mm 2 /s 3 , are 
quoted for the very deep sea. This extends the ranges to 
4 /dns ~ 0.01—50 and $dns ~ 0.03—10. 

Method. Eqs. (1,2) can be solved by iteratively refining 
approximations for the path a particle takes through the 
flow [24, 29, 30]. This results in expansions of steady- 
state averages in powers of Ku and allows to determine 
how the remaining parameters (<£>, T, and A) affect pref¬ 
erential sampling and small-scale clustering. The details 
of this calculation are given in the Supplemental Ma¬ 
terial [31]. Here we outline the essential steps. First, 
to consistently track the orders in the expansion we de- 
dimensionalise t' = t/r , r' = r/p, u' = u/u 0 . Second, we 
expand the dynamics of the vector n in powers of Ku: 

OO 

n(t') = ^ n q {t') Ku 9 . (3) 

< 2=0 

Inserting this ansatz into (1,2) and identifying terms of 
order Ku 9 yields equations for n q that can be solved in 


terms of n p for p < q. The lowest-order solution in Ku 
is just n 0 = 9 This yields a lowest-order deterministic 

approximation for the particle position at time t'\ 

T det( i ') = r o -$9? . (4) 


Third, we expand Eqs. (1,2) in terms of deviations from 
the deterministic trajectory Sr'(t') = r'{t!) — r' det (t') pre¬ 
cisely as described in Ref. [24]. In the fourth and final 
step we average over the fluid-velocity fluctuations in the 
statistical model. In the remainder of this Letter we sum¬ 
marise the results obtained in this way. 

Preferential sampling. Consider the steady-state aver¬ 
ages of the 2 -component u z of the fluid velocity and of 
its gradient, A zz , both evaluated at the particle position. 
Analytical results for these averages are derived to lowest 
order in Ku in the Supplemental Material [31], Eqs. (S15) 
and (S16). These expressions are plotted in Fig. 1. Here 
we quote only limiting results. For small 4> we have 
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What can we learn from these analytical results? 
Eqs. (5a) and (6a) show that the particles collect in 
the sinks of the transversal flow-velocity field, Tr^ A = 
—A zz < 0. This is because gyrotaxis breaks up-down 
symmetry: when T is small the particles swim essentially 
upwards (in the e z -direction), and gather in transversal 
sinks irrespective of their shape and swimming speed. 
Simulations (Fig. la) confirm the theory. 

Motivated by Kessler’s study in pipe flows [7] the au¬ 
thors of Ref. [11] concluded that spherical particles pref¬ 
erentially sample down-welling regions also in turbulence. 
This is not in contradiction with the result discussed 
above because particles may preferentially sample dif¬ 
ferent observables. In fact Eqs. (5b) and (6b) explain 
that spherical particles are biased towards down-welling 
regions (as observed in DNS [11]), in addition to sinks in 
the transversal flow. But (6b) also shows that elongated 
particles [A > d/(d + 2)\ preferentially sample up-welling 
regions for large enough <f>. This is seen in Fig. lc which 
shows (u z )oo (Eq. (S15) in the Supplemental Material), 
for Ku = 0.1 and d = 2 as a function of 4>. Also shown 
are results of statistical-model simulations, in excellent 
agreement with theory. Fig. Id shows that the same 
conclusions hold in three spatial dimensions for Ku = 1. 
Rods sample upwelling regions when 4) is larger than (ap¬ 
proximately) unity. 



3 



<j> $ 




consider the dependence of (w~)oo 011 the reorientation 
time 2$. Fig. If shows that the statistical-model result 
becomes independent of Ku for large Ku, and that it 
reproduces the DNS results fairly well, it explains the 
£8/tk -dependence of (u z )^ of the DNS results up to a 
prefactor of order unity. This factor is due to the fact that 
fully-developed turbulent velocity fluctuations in the dis¬ 
sipative range differ from those in the statistical model: 
they are not Gaussian, more persistent, and the proba¬ 
bility of straining regions to occur is higher [24]. 

Small-scale clustering. Which mechanisms cause two 
particles caught in the same flow region to actually col¬ 
lide? This is a two-particle problem, only indirectly 
related to preferential sampling. Fluctuations in the 
separations between nearby particles are determined by 
the dynamics of the particle-velocity gradients dvi/drj. 
Small-scale clustering occurs where V • v < 0. We have 
computed (V • v}^ to lowest order in Ku. The result is 
quite lengthy [Eq. (S32) in the Supplemental Material]. 
For small $ the full expression simplifies to: 

(V • v) oo r]/u 0 ~ -Ku(<F'F) 2 .B d (A) for 4> < 1, (7) 


FIG. 1: Preferential sampling, a (Trx A),*, as sampled by the 
particles. Results of simulations of statistical model: symbols. 
Eq. (5a): lines. Parameters: Ku = 0.1, 4/ = 0.1, A = 0 (o), 
= 1, A = 0 (0), 4" = 0.1, A = 1 (□) 4> = 1,A = 1 (a). 
b The same but for Ku = 1 (theory only for large d>). c 
(u z )oo as sampled by the particles for Ku = 0.1. d Same but 
for Ku = 1. e {u z )oo/v s as a function of 4? for Ku = 0.1. f 
(u z )oo/v s versus 2§/t k for A = 0 and Ku = 1 (▼), Ku = 2 
(*), Ku = 5 (◄), Ku = 10 (►). Hollow markers show DNS 
data from Fig. 3d in Ref. [11] (at ReA = 64). All data are for 
values of 4? from the small-"!? plateau observed in DNS [11], 
and also in the statistical model (for Ku = 0.1 this plateau is 
shown in panel e, for Ku = 1 in Fig. SI in the Supplemental 
Material). Panels a,c,e are for d = 2, panels b,d,f for d = 3. 


It is remarkable that the shapes of the curves at Ku = 
0.1 are very similar to those at large Ku. This means 
that the small-Ku theory qualitatively explains what is 
observed in the statistical-model simulations at large Ku 
and in DNS for spherical particles [11]. 

In the limit of large 4? particles swim rapidly upwards 
and experience the flow as a white-noise signal (as in 
rapid gravitational settling [30]). This limit is universal, 
particles in any homogeneous, isotropic and incompress¬ 
ible flow show preferential sampling according to Eq. (6). 
This means that the small-Ku theory should describe re¬ 
sults of statistical-model simulations at Ku = 1 quanti¬ 
tatively for large 4?. This is confirmed by Figs. lb,d. 

For small 4? DNS [11] show that the average of u z is 
proportional to 4? for small 4>, so that (m 2 ) 00 /^ ) is con¬ 
stant. Eq. (5b) shows this behaviour, in good agreement 
with simulations (Fig. le). 

We conclude with a quantitative comparison of 
statistical-model and DNS results [11]. As an example 


with B d { A) = [(d+2)(d+4)-2d(d+4)A + (4+2d+d 2 )A 2 ]/d. 
Since B d ( A) > 0, Eq. (7) implies small-scale clustering. 
For spherical particles (A = 0) the quadratic dependence 
of (V • u)oo on 4>4/ was derived in Ref. [11] (and also in 
Ref. [32]): expanding Eqs. (1,2) for S3 <C r gives 

V • v ~ u s ^[-(l+A) d1u z + (1-A)(9 2 M 2 -Ait 2 )] . (8) 

Substituting A = 0 yields Eq. (6) of Ref. [11], and aver¬ 
aging Eq. (8) along particle paths results in Eq. (7). The 
factor v s 2d in (8) corresponds to one factor of 4? 4/ in (7). 
The second factor of 4> T comes from averaging the ve¬ 
locity derivatives in Eq. (8). We note that Trx A does 
not figure in Eq. (8): preferential sampling of sinks in 
the flow-velocity field perpendicular to gravity does not 
contribute to small-scale clustering, showing that the two 
effects are distinct [24]. 

Expanding the full result (S32) for large 4? gives: 

(V ■ v) oo r]/u 0 ~ - Ku 4?'F 2 E d (A) for 4? » 1, (9) 

withKd(A) = ^/tt/ 2(d+l)(d+3)(A — l) 2 /d. For spherical 
particles the 4”F 2 -dependence was derived in Ref. [32]. 

Let us now analyse the shape dependence of Eq. (7). 
The A-dependence of B d ( A) explains that rods (A = 1) 
cluster less than spheres (A = 0), consistent with the 
DNS results reported in Ref. [22]. But when gyrotaxis 
is weak, spheres are essentially randomly oriented, unlike 
neighbouring rods that are aligned by turbulent shears. 
In this limit motile rods must cluster more than spheres. 
This is demonstrated below, but it is not captured by 
Eqs. (7) and (9) which must fail for large 4> because the 
limit T —> oo is singular, and due to the occurrence of 
singularities in the dynamics of the gradients of n at large 
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but finite values of (Supplemental Material). The first 
caveat also applies to Eqs. (5) and (6). 

Fractal dimension. DNS show [11] that the small-scale 
spatial patterns of motile gyrotactic organisms are frac¬ 
tal. This may substantially enhance their encounter rates 
[33]. We analyse the fractal patterns for finite but small 
Ku, in two dimensions. We expect qualitatively the same 
result in three dimensions. The fractal patterns are char¬ 
acterised by ‘Lyapunov exponents’ Ai and A 2 

Ai = lim f -1 In —— and A 1 + A 2 = lim f -1 In . (10) 

t —>00 _R(o) t —>00 -4(o) 

These exponents quantify the expansion (contraction) 
rates of the distance R(t) between two initially nearby 
particles, and of the area element A(t) spanned by the 
separation vectors between three nearby particles. The 
fractal Lyapunov dimension is defined by [24, 34] 

dh = 1 — A 1 /A 2 , (11) 

assuming Ai > 0 and Ai + A 2 < 0. When d\, < 2 fractal 
clustering occurs. To evaluate d\, we use Ai + A 2 = (V • 
u)oo, Eq. (S32), and compute Ai to order Ku 4 . The result 
is lengthy, in the Supplemental Material [31] we quote 
the result to order Ku 2 , Eq. (S31). To this order it is 
independent of T and A. For small values of <f> Eq. (S31) 
simplifies to Air ~ Ku 2 (l — 3 4> 2 ) for d = 2. Together 
with (7) this implies Al = d — c?l ~ <f ,2 4' 2 consistent 
with the results of Refs. [11, 32] for spherical particles. 
Fig. 2a shows the analytical result for d^ as a function of 
<f>. It is in good agreement with numerical simulations of 
the statistical model (d = 2) for Ku = 0.1. We see that 
spherical organisms cluster more than rods. As explained 
above this is expected for strong gyrotaxis. 

But when the effect of the gravitational torque is small 
then prolate organisms cluster more: in the absence of gy¬ 
rotaxis, rotational symmetry ensures that active spheri¬ 
cal particles remain uniformly distributed, but rod-like 
particles show fractal clustering. Panel b in Fig. 2 
demonstrates this cross-over. It shows dh for 4> = 1, 
Ku = 1 as a function of 4/. We arrive at qualitatively 
similar conclusions by numerically computing the frac¬ 
tal correlation dimension c? 2 . But the numerical values 
found for differ from This shows that the spatial 
distribution is multifractal [24]. 

Conclusions. First, our statistical-model calculations 
explain how the dynamics of gyrotactic motile micro¬ 
organisms depends on the dimensionless parameters of 
the problem: A (shape), <f> (swimming speed), and T 
(reorientation time). Second, we find that the particles 
tend to preferentially sample positive values of A zz , cor¬ 
responding to sinks in the transversal flow, regardless of 
shape. We predict that this must also be observed in 
DNS, it is simply a consequence of the fact that gravity 
breaks the symmetry of the problem. At the same time 
our calculations show that spherical particles are more 



FIG. 2: a Fractal dimension deficit A l = d — for d = 2, 
Ku = 0.1. Numerical simulations of statistical model for 4> = 
0.1, A = 0 (©), 4' = 1, A = 0 (0), 4- = 0.1, A = 1 (□), and 4< = 
1,A = 1 (a). Theory [Eqs. (11), (S31), and (S32)] including 
the Ku 4 -contribution to Eq. (S31): solid lines. Asymptote oc 
(<j>4>) 2 : dashed lines, b Numerical simulations, c(l for d = 2, 
$ = l, Ku = 1, A = 0 (o), 0.2 (□), 0.4 (0), 0.6 (a), 0.8 (▼), 1 
(★). 

often found in regions where u z is negative, explaining 
the behaviour found in DNS [11]. But our calculations 
also predict that rod-like particles preferentially sample 
up-welling regions of homogeneous isotropic flows such as 
turbulence, provided that they swim fast enough. Third, 
we have analytically computed how the degree of small- 
scale spatial clustering depends on particle shape. This 
is important because small-scale fractal clustering may 
enhance particle-encounter rates. We find a transition 
that we predict to be observable in DNS as well: when 
gyrotaxis is strong (small T) oblong particles cluster less 
than spherical ones, while at large T the opposite is true. 

Our calculations also show that singularities in the mo¬ 
tion of nearby micro-organisms occur, much like ‘caus¬ 
tics’ for heavy particles in turbulence [35-38]. We pre¬ 
dict that such singularities must also be observed in the 
DNS of gyrotactic microswimmers in turbulence. It is 
of interest to estimate how often the singularities occur 
because their effect may modify the predictions of phe¬ 
nomenological models for encounter rates [39]. 

The analytical results obtained in this Letter were de¬ 
rived for small Ku (or large <f>). But we have shown 
that our analytical results and the corresponding mecha¬ 
nisms qualitatively explain what is observed in DNS, and 
explain also the results of statistical-model simulations 
at large values of Ku. We find fairly good quantitative 
agreement between our statistical-model calculations and 
DNS results for fully developed turbulence. To achieve 
even better quantitative agreement with the DNS would 
require to account for the universal non-Gaussian small- 
scale fluctuations of fully developed turbulence [26]. 

But the fluctuations of the unsteady ocean are nei¬ 
ther fully-developed turbulent, nor are they Gaussian. 
Therefore the fact that the much simpler Gaussian sta¬ 
tistical model explains the dynamics observed in DNS 
of fully developed turbulence [11] shows that the ana¬ 
lytical theory (and the underlying mechanisms) describe 
robust behaviour, that must be taken into account in 
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the analysis of patchiness and encounter rates of motile 
micro-organisms in the ocean. 
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In this Supplemental Material we outline the steps that are needed to derive Eqs. (5,6) in the Letter [SI]. This 
Supplemental Material also describes how we computed an analytical expression for the Lyapunov dimension shown 
in Fig. 2a of the Letter, and it contains a Supplemental Figure (Fig. SI). 


PREFERENTIAL SAMPLING 


The results described in the Letter [SI] rest on the perturbation method described in Refs. [S2-S4]. In this Section 
we summarise how we have applied this method to the equations of motion (1,2) in the Letter, in order to compute 
the average of the fluid-velocity held and its gradients along particle paths. These averages determine to which extent 
the particle preferentially sample certain huid-velocity configurations. 

In the dimensionless variables t' = t/r,r' = r/p,u' = u/uq,uj' = ujt the equations of motion read: 

dr' 

— — = Ku u + dm (SI) 

dr 

^ =U}' An (S2) 

u/ = — (n A g)/{ 2H/) + Ku + Ku An A [§ / (r / (t / ), t!)n\ . (S3) 


The dimensionless parameters Ku = uqt/t /, T = 38 and $ = v s t/t] are chosen so that all terms containing the 
fluid velocity u or its gradients are proportional to Ku. This implies as we shall see that a perturbation expansion 
about the deterministic (u = 0) solution is equivalent to a perturbation expansion in the Kubo number Ku. 

To expand Eqs. (SI) and (S2) in the Kubo number we first require a Ku-expansion of the orientation vector n: 


n(t') = ^2n q (t') Ku 9 , 

9=0 


(S4) 


where n q (t') are time-dependent expansion coefficients. Inserting this ansatz into Eq. (S2) and identifying terms of 
order Ku 9 we find that the expansion coefficients must satisfy: 

Sl) n j(t') - 6 qfi gj + M'(r , {t , ),t , )n q _ 1 {t') 

3=0 

“X X ( n i(0 ■M'(r'{t'),t')n k {t'))n q ^ j - k - 1 {t') . (S5) 

j=0 k —0 


Here B' = O' + AS', and S' and O' are the symmetric and antisymmetric parts of the matrix A' of fluid-velocity 
gradients. Eq. (S5) can be solved for n q in terms of integrals over solutions n p with p < q. The lowest-order solution 
in Ku, n 0 (t'), is deterministic. It does not depend on B, and is thus independent of the stochastic fluid-velocity 
gradients. This deterministic solution no(t / ) depends on the initial orientation in a complicated way. But in the limit 
of large times the dependence on the initial orientation must drop out (provided that 4/ 0). Using this fact we 

obtain no(t' —> oo) = — g. We have verified explicitly that the dependence on the initial orientation drops out for all 
averages evaluated in the Letter, by asymptotic expansion of n 0 (t') for large values of t and by numerical evaluation 
of the exact solution for n 0 (t'). 

Using no(t') = — g we find the following solution for the expansion coefficients n q (t'): 


n„ 


(f) = s q , 0 g + f d U-*')/^)/^) + (e^-*')/* - e (*i-*')/(2*)) (/ ( 4 / 1 ). g)g 

Jo 1 


(S6) 
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with 


9-1 


9-19-J-l 


fit') = ^ (' n q-o(t') ' g)nj{t') + W {r' -J2Y1 ' B V(0>0 n fc(0)™9-i-fc-i(0 • 


j=i 


j —0 /c—0 


(S7) 


Now we insert the resulting solution for n(t') into the equation of motion (SI). This yields: 

dr-' 00 

— =Ku«'(r'(t')d') «,,(£') Ku 9 . (S8) 

9—0 

This equation has the implicit solution 

pt' oo 

r \t') = r' det {t') +Ku / dt , 1 [u'(r , (t , 1 ),t , 1 )+$^n (? (t , 1 )Ku 9 " 1 ] . (S9) 

J 0 g=l 

Here r f ( et (t') denotes the deterministic solution, obtained by letting Ku = 0 in Eq. (S8) and using n 0 {t') = —g: 

»det(0 = r \0) - Ku A>gt '. (S10) 

Now consider an actual particle path r'(t'). The difference 5r'{t') = r'(t') — rj et (t / ) between this path and its 
deterministic approximation is obtained from Eq. (S9). We see that 5r'(t.') ~ Ku. This ensures that 8r'(t') is small if 
Ku is small enough, and it allows us to expand the fluid velocity u'(r'(t'),t') in Sr' around the deterministic trajectory 

»det (*)’■ 

u'(r'(t'),t') = u' (r^ et (<:'), t') + [Sr'(t') ■ V'] v! {r' Aet {t'), t') + ... . (Sll) 

Inserting this expansion and the corresponding expansion for B '{r'{t!),t') in (S7) into Eq. (S9) gives an expression 
for Sr'(t') in terms of t') and its derivatives, analogous to the expressions discussed for the inertial-particle 

problem in Ref. [S4]. 

Finally, to evaluate the average (u)^ we iteratively substitute Sr'(t') into Eq. (Sll), and cut off the resulting series 
at the desired order in Ku. To first order in Ku we find for example: 


= u'{r' det (t'),t') + KuJ dt , 1 |[w , (r , (t , 1 ),ti) +$n 1 (t , 1 )] • V'ju'(rd et (t , ),t'), (S12) 


with 


n i 


(t'i) = f d t' 2 e {t ' 2 t ' l)/(2 ' I ' ) { -M'(r'{t' 2 ),t' 2 )g+ [g-M’{r’(t' 2 ),t' 2 )g]g} . 


(S13) 


Up to this point the expressions are valid for fluid-velocity fields with arbitrary distribution. Now we make use of the 
fact that the fluid-velocity field is assumed to be Gaussian in the statistical model with, zero mean, Gaussian spatial 
correlations, and with exponential time correlations [S4]: 


(u'(r i,tiK-(r',4)) = 


1 


d(d- 1) 


[Said ~ 1 - ( r i - r' 2 ) 2 ) + (ri - r' 2 Ur\ - r' 2 ) 3 ] e -W^) 2 /2-|t' 1 -4l . 


(S14) 


The form of the prefactor follows from isotropy, homogeneity, incompressibility, and normalisation (tt'(0, 0) 2 ) = 1. 
Taking the average of Eqs. (S12) and (S13) we find: 
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Ay/2 $ 3 4 
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(S15) 


with C x = A + 1 + (d + 1)(A — l)x 2 and T[x\ = ypKe 3 ? erfc[x]. In the same way we obtain the average of the 
fluid-velocity gradient matrix A': 


(A')oo 


Ku$ Sjj - dg-jgj f 4T + 1 

d 1 - d l° 0 4$ 3 T 


'hC<j> ^ 

V^ 4 


r 1 i 
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8\/2 $ 4 A! 2 


[■ 2’k + 11 
-V / 8$’K 


}. (S16) 
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The ^-components of Eqs. (S15) and (S16) are plotted in Fig. 1 in the Letter [SI]. The asymptotic behaviour for 
small values of $ is obtained by series expansion in <f>: 


«>c 


Ku $' 


- Ku$ 


d(l - A) + 2(A + 2) >F(4>F + 1) 
d (2<F + 1) 2 ’ 

d{ l-A) + 2 ^ 


2T + 1 


(517) 

(518) 


These equations correspond to Eqs. (5a) and (5b) in the Letter. 

A corresponding series expansion for large $ of Eqs. (S15) and (S16) gives 


(A'zz) c 

«>= 


Ku d +1 


(1-A h/j, 


<F 2d 
Ku (d (A - 1) + 2A) 


$ 


2d 


(519) 

(520) 


These are Eqs. (6a) and (6b) in the Letter. We note that this result can also be found by an expansion for large $ 
of the original equations for an arbitrary homogeneous, isotropic and incompressible flow, for arbitrary values of Ku. 
Because this limit is universal it is not necessary to assume that Ku is small. 
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SMALL-SCALE CLUSTERING 


The small-scale spatial patterns of gyrotactic micro-swimmers are fractal [S5]. Fractal small-scale clustering is 
quantified by fractal dimensions (see Ref. [S4] for a review). There are different fractal dimensions, characterising 
different geometrical aspects of the spatial distribution. In the Letter we show results for the fractal Lyapunov 
dimension in d = 2 spatial dimensions. Eq. (11) in the Letter shows that this dimension is defined in terms of the 
Lyapunov exponents: 

A ‘ TS A” H “ d Oi+A 2 )r Si lim (S21) 


The exponent Ai gives the rate which the separation R'(t') between two initially nearby micro-swimmers decreases 
or increases exponentially at long times. The sum Ai + A 2 determines the rate at which the infinitesimal area A!(j?) 
of the parallelogram spanned by the separation vectors between three initially nearby micro-swimmers decreases or 
increases at long times. The problem of computing small-scale fractal clustering in d dimensions is a d + 1-particle 
problem. 

We solve this problem using the method described in Ref. [S4]. Linearising the equation of motion Eq. (SI) we see 
that the separation between two initially nearby micro-swimmers evolves according to 


d R! 
df 


Ku (R ■ UR) R ', 


(S22) 


where R denotes the unit vector along the separation direction between the two micro-swimmers, and U = A' + $ ¥' 
is the particle-velocity gradient matrix. Here ¥' is the matrix of orientation gradients with elements = dn t /dr' :j . 
It follows from Eq. (S22) that the exponent Ai is given by 


Air = Ku (R ■ UR) 00 . (S23) 

In two spatial dimensions the sum of Ai and A 2 is simply 

(Ai+A 2 )r = Ku(TrZ , )oo- (S24) 


To evaluate these averages we require the time evolution of R and ¥'. This follows from linearisation of Eqs. (SI) 
and (S2): 


^ = Ku [Z'jRj - (RjZ'MRi ], (S25) 

d Y' l 

~^T = ^\3k Y kj n i + n k g k Y[.\ + Ku[(djB' ik )n k + B\ k Y^ - 2A(n k S kl Y/ j )n i - A {n k (djS' kl )m)ni 

- A (n k S' kl niW j - Y' k (A' kj + <t>YQ] , (S26) 

where repeated indices are summed over. We solve Eq. (S26) in terms of a series expansion 


Y'(t') =£¥;(*') Ku* . 

q=0 

Inserting this expression into Eq. (S26) we find the following expression for the expansion coefficients: 

rt' 


(S27) 




dt( 


e (ti-t')/(2®) i r. .ft) + ( e «-*')/^ _ eK- t 'VW)g i g k F kj (l' 1 ) 


(S28) 


where we the initial condition was put to zero because the corresponding terms do not contribute to steady-state 
averages. The coefficients F;,,j(t') are given by 


9 -i 


Fiji 1 ') = 'Yl[h Y r-kj n q-r-,i + n q - r -k9k Y r-,ij + S qA( d jB'i k (r'(t'),t'))n k - A n k {djS' kl (r'(t'),t')ni)ni 


r—0 


+ B' ik {r\t'), t!)Y^ - 2Kn k S' kl {r\t'),t')Y{ j n i - A n k S' kl (r'{t'), t'^Y^ - Y{ k A' kj (r'- $Y{ k Y kj ]. (S29) 
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Eq. (S25) can be solved implicitly as 


Mf) = Ku f d A{z'ii{A)RAA) - [RAW k {A)h{A)]iu{A)} ■ 


(S30) 


By substituting R and ¥' in the right-hand side of Eqs. (S29) and (S30), Eqs. (S29) and (S30) can be recursively 
iterated to any desired order in Ku. Expanding the fluid velocity and gradients of the fluid velocity as in the previous 
Section gives equations for R and ¥' in terms of the fluid velocity evaluated along the deterministic trajectories. 
Averaging gives the Lyapunov exponents to order Ku 2 


Air = 


Ku 2 


: { [(d-3) c f >2 +{3 —d+(15 —4d+d 2 ) <^ 2 ){(R • g) 2 )-2(l + 5 <S> 2 )((R ■ g) 4 


$+4?f(3-d)$ 2 (S31) 


d(d- 1)$ 5 IL V “ - - / -V- — /j ^ L V 

(3-2 d+d 2 ) 4> 4 —(3 —d+(18-5d—d 2 ) $ 2 +(9-2 d+d 2 ) § 4 )({R ■ g ) 2 ) + (2 + 12 $ 2 +6 $ 4 )((J k ■ g) 4 ) E[1/{V2 $)]} , 

{(A + 1) [C$ + $ 2 (d(A - 1) - 2)] $-[C| + 2(d+l)(A-l) 2 $ 4 ] J'[l/(v / 2$)]/v / 2}. (S32) 


(Ai + A2 )t — 


Ku 2 ^ 2 


d$ d 


The steady-state moments {(R ■g) 2p ) oc in Eq. (S31) characterise the anisotropy of the spatial patterns. The moments 
are given by Eq. (8) derived in Ref. [S3] for particle pairs settling in a turbulent aerosol, upon replacing G in that 
equation by — <!>. The fact that the moments agree is a coincidence, only valid to lowest order in Ku. While 
the problems are superficially related (to lowest order particles move on a deterministic trajectory r^ et (t l ) through 
turbulence) the details are quite different. Using the above expressions for the Lyapunov exponents (with Ai extended 
to order Ku 4 ) we obtain the curves plotted in Fig. 2a in the Letter. 

Series expansions of Eq. (S32) for small and large values of <f> give Eqs. (7) and (9) in the Letter 

(V • v) oo ri/u 0 ~ - Ku ($4') 2 R d (A) for $ < 1, (S33) 

(V • i^ooTy/wo ~ — Ku $4' 2 K ( i(A) for 4> 1. (S34) 


The shape-dependent prefactors Bd( A) and Ed{ A) are found to be: 

B d ( A) = [(d + 2)(d + 4) - 2d(d + 4)A + (4 + 2d + d 2 )A 2 ]/d (S35) 

E d (A) = V / V2(d+l)(d + 3)(A-l) 2 /d. (S36) 


The corresponding limits for Ai are more subtle due to the complicated dependence of the moments ((R ■ g ) 2p )oo on 
4>. For small values of $ we have ((R ■ g ) 2p )oo ~ (2 p — l)!!/(2 p p!) [S3] and consequently 


Air 


d 2 + d — 3 
d(d — 1) 


Ku 2 $ 2 . 


(S37) 


This equation is quoted in the Letter for d — 2. 

The results derived here are obtained for general values of $, tp, and A. It is important to note that the limit of 
large is singular. In this limit the dependence on the initial orientation n(t = 0) does not decay fast enough for it 
to be disregarded. This explains, as pointed out in the Letter, that the theory must fail for very large T. 

We also note that Eqs. (S31) and (S32) are derived assuming that the elements of ¥' are small enough, allowing us to 
expand Eq. (S26) in Ku. This may fail when singularities occur in the dynamics of Y'. Elements of ¥ may repeatedly 
diverge to —oo. These singularities are analogous to caustics in inertial-particle dynamics [S4] and correspond to 
instances where the n-field becomes multi-valued so that nearby particles can have very different orientations. The 
perturbation expansion leading to Eqs. (S31) and (S32) is expected to diverge when caustics are frequent, precisely 
as in the inertial-particle problem (this point is discussed in detail in Ref. [S4]). 
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SUPPLEMENTAL FIGURE 


Fig. SI shows that (u z )oo/v s approaches a constant as $ —> 0. The parameters are Ku = 1, d = 3, T = 1, and 
A = 0. The Figure demonstrates that the plateau forms also for larger Kubo numbers, not only for small Kubo 
numbers (Fig. le in the Letter shows a corresponding plateau for Ku = 0.1). 



FIG. SI: Shows (u z )oo as a function of $ for Ku = 1, d — 3, 'F = 1, and A = 0. The data used in this Figure are identical to 
the data shown in Fig. Id in the Letter. The vertical dashed line shows the value of <1> used for the Ku= 1-data (▼) in Fig. If 
in the Letter. 
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